Self- Assembly of Monatomic Complex Crystals and Quasicrystals 
with a Double- Well Interaction Potential 
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For the study of crystal formation and dynamics we introduce a simple two-dimensional 
monatomic model system with a parametrized interaction potential. We find in molecular dy- 
namics simulations that a surprising variety of crystals, a decagonal and a dodecagonal quasicrystal 
are self-assembled. In the case of the quasicrystals the particles reorder by phason flips at elevated 
temperatures. During annealing the entropically stabilized decagonal quasicrystal undergoes a re- 
versible phase transition at 65% of the melting temperature into an approximant, which is monitored 
by the rotation of the de Bruijn surface in hyperspace. 

PACS numbers: 61.50.Ah, 02.70.Ns, 61.44.Br, 64.70.Rh. 



Self-assembly is the formation of complex patterns out 
of simple constituents without external interference. It 
is a truly universal phenomenon, fundamental to all sci- 
ences Although usually the constituents interact 
only locally, the result is well-ordered over long distances, 
sometimes with a high global symmetry. In the process 
of crystallization, particles (atoms, molecules, colloids, 
etc.) arrange themselves to form periodic or quasiperi- 
odic structures. Here we are interested in structurally 
complex phases. Examples are metallic crystals with 
large unit cells hundreds or thousands of atoms - known 
as complex metallic alloys @] ■ Some consequences of the 
complexity are the existence of an inherent disorder and 
the formation of well-defined atomic clusters Re- 
lated alloys differ by the cluster arrangement. In the 
limit of infinitely large unit cells non-periodic order like 
in quasicrystals [4j is obtained. However self-assembly of 
complex phases is not unique to alloys. Recently micellar 
phases of dendrimers were observed to form a mesoscopic 
quasicrystal [5| , and there are indications that quasicrys- 
tals exist in monodisperse colloidal (macroscopic) sys- 
tems [fjj]. Since the interaction between colloidal par- 
ticles can be tuned in various ways, these systems are 
well-suited for experiments investigating self-assembly in 
dependence of the potential shape. 

All of the previous examples have in common that the 
crystal growth can be modeled with effective pair poten- 
tials, which is a prerequisite for simulating self-assembly 
of a large number of particles on a computer. The first 
such simulations were conducted by Lancon and Billard 
in two Q and by Dzugutov in three dimensions [||. In 
the latter work the system was chosen monatomic to fa- 
cilitate computation and separate chemical from topo- 
logical ordering. It is well-known that many common 
pair-potentials like the Lennard-Jones (LJ) potential fa- 
vor close-packed ground states. To force the formation of 
alternative structures the Dzugutov potential is endowed 
apart from its LJ minimum with an additional maximum. 
Although the potential was originally tailored to lead to 
a glassy state [8( , it stabilizes the er-phase at low temper- 



ature @. Later, similar potentials were used to demon- 
strate the formation of a dodecagonal quasicrystal in two 
dimensions 10], which on closer inspection was identified 
as a periodic approximant. In fact it is often difficult 
to distinguish quasicrystals and closely related periodic 
complex crystals due to their structural similarity. 

The relation between an interaction potential and 
its energetic ground state is a fundamental problem of 
physics. It can be approached by two methods: The 
direct method starts from a given parametrized set of 
potentials and studies the resulting structures as a func- 
tion of the parameters (and temperature/pressure). An 
example is the hard core plus linear ramp model with the 
ramp slope as single parameter ll( . The inverse method 
tries to find an appropriate potential that stabilizes a 
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given structure via optimization 
cently by Rechtsman et al. to find potentials for various 
lattices [13(. In this Letter we apply the direct method 
to a potential of the form 
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which we denote Lennard-Jones-Gauss (LJG) potential. 
For most values of the parameters V(r) is a double- well- 
potential with the second well at position ro, of depth e 
and width a (Fig. [T]) . We note that the general form of 
pair potentials in metals consists of a strongly repulsive 



core plus a decaying oscillatory (Friedel) term [14J. A 
LJG-potential can be understood as such an oscillatory 
potential, cut off after the second minimum. 

In the following we restrict the parameter space by 
fixing er 2 = 0.02. The T — phase diagram in the ro-e- 
plane is determined in two steps: First, candidate ground 
state structures are obtained from annealing simulations. 
Next, a defect free sample of each candidate structure 
is relaxed with a conjugate gradient algorithm. During 
the relaxation particle movements and adjustments of the 
simulation box dimensions are allowed. The stable struc- 
tures (within the candidates) are the ones with the low- 
est potential energies. Simulations were carried out by 
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FIG. 1: (color online) LJG-potential for e = 1.1, a 2 = 0.02, 
and r = 1.3 (Hex2 in Fig.©, 1.4 (Sqa), 1.5 (Pen), 1.6 (Hexl). 



solving the equations of motion with molecular dynam- 
ics. Periodic boundary conditions and a Nose-Hoover 
thcrmo-/barostat for constant temperature T and con- 
stant pressure P = were used. A cut-off was applied to 
the LJG-potential at r = 2.5. 

We performed 5000 annealing simulations with a sam- 
ple of 1024 particles using parameters located on a fine 
grid: r e [1.11,2.10], Ar = 0.01 and e e [0.1,5.0], 
Ae = 0.1. During each run the temperature was lowered 
linearly over 2 • 10 6 molecular dynamics steps starting 
from the liquid state. A typical final particle configura- 
tion consists of several well-ordered grains whose struc- 
ture is analyzed in direct and in reciprocal space. The 
phase observation regions are displayed in Fig. [21 We 
found hexagonal phases with nearest-neighbor distance 
close to 1.0 (Hexl) and close to r (Hex2), a square lat- 
tice (Sqa), a decagonal (Dec) and a dodecagonal (Dod) 
quasiperiodic random tiling (RT) , a phase built from pen- 
tagons and hexagons (Pen), a rhombic lattice (Rho), and 
finally a honeycomb lattice (Hon). It can be seen in Fig. [2] 
that the two hexagonal phases are connected around the 
phase Sqa. Across the line between C and the phase Sqa 
there is a rapid increase in the hexagonal lattice constant. 

The phases resulting from the annealing simulations 
are possible energy ground states. Further possibilities 
are approximants, which we constructed from the tiles in 
Fig. [21 The choice of approximants is restricted by the 
following observations at low temperature: (i) The D- 
tiles have lowest energy and their number is maximized 
(see below), (ii) The Sh-tiles and two neighboring Sq- 
tiles are avoided. Together, the approximants and the 
phases from annealing simulations were used as initial 
structures for numerical relaxation. In the phase dia- 
gram (Fig. [3]) five complex crystals are stable: the phase 
Pen, the decagonal approximant (Xi), which is a periodic 
stacking involving D-tiles, and three dodecagonal approx- 
imants: the cr-phase (Sigl) and two modifications (Sig2), 
(Sig3). Here we use the term "complex" since the lat- 
tice constants are larger than the potential cut-off, which 
means that the unit cells have to be stabilized indirectly 
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FIG. 2: (color online) Phase observation regions after anneal- 
ing simulations. Tilings of the unit cells are shown. The qua- 
sicrystal tiles are: (P)entagon, (H)exagon, (N)onagon, (U)- 
tile, (D)ecagon, (Sq)are, (Tr)iangle, (Sh)ield. 
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FIG. 3: (color online) Phase diagram of the LJG-potential at 
T = with a 2 = 0.02. Four approximants have been found. 
Phason flips correspond to local changes in the tilings. 



by geometric constraints. Quasicrystals are not energet- 
ically stable at T = 0. The phase boundaries are slightly 
displaced compared to those from annealing simulations 
due to metastability. 

The location of the stability regions can be understood 
from the near-neighbor configuration. Local 71-fold order 
is stabilized for ro ~ 2cos(7r/n). Fig. [3] confirms this ex- 
cept for local five- fold order, that is found at ro ~ 1.47 as 
a result of the competition with the close-packed phase 
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Pen 


0.8981 


a = 2.62, b = 2.50 
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0.7617 


a = 4.24, b = 4.24 
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1.0718 
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1.0829 


a = 4.73, b = 2.73 
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TABLE I: Structure details of the ideal tilings of the complex 
crystals. For comparison the quasicrystals are included. 



Hexl. Structure details of the complex phases are col- 
lected in Tab. [TJ The phases Xi and Dec have a sur- 
prisingly low density, which is only 65% of the density 
of Hexl. They are locally very similar, their potential 
energies differ by less than 0.7%. 

We have studied the decagonal quasicrystal in longer 
simulations using the parameters r = 1.52, e = 1.8, and 
a 2 = 0.02 (indicated by a cross in Fig. [5]). With these 
parameters, the decagonal quasicrystal is assembled with 
few defects at elevated temperatures. A simulation of a 
large sample, 10000 particles, was initiated in a random 
configuration at T = 0.50, which is close to the melt- 
ing point Tm = 0.56 ± 0.02. In the following, the peri- 
odic boundary conditions are turned off to allow phason 
strain relaxation. At the beginning the particles quickly 
formed an amorphous state with local decagonal order. 
After about 10 5 molecular dynamics steps multiple grains 
with the quasicrystal started to grow. The bigger ones 
increased their size until at ca. 10 7 steps only one single 
grain remained. We continued the simulation up to 5-10 7 
steps, healing out point defects (vacancies, interstitials) 
and improving the quasiperiodicity. At the end the sam- 
ple was quenched to T = and relaxed. The diffraction 
image (Fig. [4| shows a perfect decagonal symmetry with 
long-range order. There is a weak pattern of intrinsic 
diffuse scattering due to the randomness. 

In the simulation the dynamics is dominated by parti- 
cle jumps over the short distance Ar = 0.6, called phason 
flips, which transform energetically comparable configu- 
rations into another (see Fig. [3]). Each such configuration 
can be mapped to a tiling and embedded as a discrete de 
Bruijn surface in a five-dimensional hypercubic lattice 
by noting that the tiling vertices are integer multiples of 
the five basis vectors e n = (cos(27rn/5), sin(27rn/5)). Al- 
though the average orientation of the surface is fixed by 
the decagonal symmetry, phason flips lead to local fluc- 
tuations h (r) in "perpendicular space" resulting in a 
phason strain Xij — Vj/d - . The ensemble of all accessi- 
ble configurations forms a cntropically stabilized random 
tiling with a phason elastic free energy density of the 
general form f(T, X ) = Xi(T) Xl 2 + A 2 (T)x 2 2 , where X i, 
X-2 are quadratic forms in Xij, an d ^2 independent 



FIG. 4: (color online) Diffraction image of the large 10000 
particle sample with the decagonal quasicrystal. 

phason elastic constants flEl ]. 

According to the T — phase diagram a transforma- 
tion to the phase Xi can occur during annealing [l6j |. 
This is achieved by (i) a collective rearrangement of the 
tiling induced by a global change of the de Bruijn sur- 
face orientation and (ii) a damping of the local phason 
fluctuations. For (i) to happen a huge number of phason 
flips is necessary, which makes the transition extremely 
slow. We performed long simulations over 10 9 steps with 
1600 particles. At intermediate temperatures, T < 0.40, 
a reversible change in the tiling was found: The density 
of D-tiles increased and they arranged preferably close- 
packed in rhombs, characteristic for the approximant Xi. 
However as shown in Fig. some defects and stacking 
faults were still present in the rhomb super-tiling. 

At low temperatures, T < 0.30, the flip frequency was 
too slow to reach equilibrium with molecular dynamics. 
Hence we turned to a Monte Carlo algorithm, which al- 
lowed sampling the full temperature range from T = 0.5 
down to zero and back up. As elementary step a random 
displacement inside a circle of radius 0.7 was used. The 
large radius allows both local relaxation and phason flips. 
The squared average phason strains Xi> X2 are indica- 
tors of the de Bruijn surface orientation and thus order 
parameters for the phase transition. The results in Fig. [6] 
show a reversible transition at T c = 0.37 ± 0.03. Above 
T c there is a strong increase in the density of D-tiles, 
which then slows down below T c . The phason strains 
X\ and xi fluctuate and switch from zero average at 
T > T c to finite values including a small hysteresis. We 
note that with periodic boundary conditions or in large 
samples phason flips alone cannot change the decagonal 
symmetry efficiently. 

In the case of the Sig-phases, phason flips are not pos- 
sible. They only occur in combination with Sh-tiles as 
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FIG. 5: (color online) Particle configuration in molecular dy- 
namics simulation at T = 0.30. The D-tiles are arranged 
preferably in a rhomb super-tiling. 



stability. In contrast, the flip frequency of earlier mod- 
els @, HI is much lower. Another example, a repulsive 
term plus two negative Gaussians has a qualitatively sim- 
ilar phase diagram, although additional phases appear. 
Further details will be presented elsewhere. 

In conclusion, we have shown that systems with two 
competing nearest-neighbor distances can have a much 
more complicated phase behavior than what is known 
for single minimum potentials. Quasicrystals and com- 
plex crystals appear naturally in such systems as an at- 
tempt to maximize local particle configurations with non- 
crystallographic symmetry. In the presence of phason 
flips, entropic contributions to the free energy play an 
important role in the thermodynamic stability. 
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FIG. 6: (color online) Monte Carlo simulation of the phase 
transition between the decagonal RT and the approximant Xi. 
The data has been averaged over intermediate temperature 
intervals. Open symbols: cooling; full symbols: heating. 

depicted in Fig. [3l Even though Sh-tiles are not seen 
in the ground states, they are present in equilibrium at 
higher temperature as structural defects. The flip mech- 
anism then differs from the one in the decagonal RT [13] • 
Finally we comment on the value of the parameter a 
and on other double-well potentials. The dynamics of 
the complex phases is controlled by the potential hill be- 
tween the minima. A high potential hill leads to a low 
phason flip frequency and slow phase transitions, at least 
in molecular dynamics. On the other hand, a too low po- 
tential hill does not stabilize complex phases. The phase 
behavior is quite robust against small changes in the po- 
tential. Phase diagrams resembling Fig. [3] have been ob- 
tained for different values of a. Choosing a 1 = 0.02 con- 
stitutes a compromise between high flip frequency and 
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